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Abstract 

Internal kink instabilities have been studied in straight tokamak geometry employing an elec- 
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r*j \ tromagnetic gyrokinetic particle-in-cell (PIC) code. The ideal-MHD internal kink mode and the 



Of 



04 



collisionless m = 1 tearing mode have been successfully simulated with the PIC code. Diamagnetic 
effects on the internal kink modes have also been investigated. 
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I. INTRODUCTION 



The internal kink instability is a fluid fixed-boundary mode which can, however, be 
strongly affected by kinetic effects. It occurs if the safety factor q(r c ) = 1 at some position 
r c inside the plasma column jj-y]. The fluid internal kink mode is always unstable in screw 
pinch geometry, but it can be partially stabilized by toroidal corrections (4]. Nevertheless, 
it plays an important role in tokamak operation (in particular during "sawtooth" 5| and 
"fishbone" Q| activity). In the fluid picture, the internal kink mode is destabilized by the 
gradients of the parallel current and of the plasma pressure. Physically, the kink instability 
is a tilt and a shift of the plasma column. It can be unstable because the field line bending 
vanishes at the position r c and fails to compensate the magnetic force associated with the 
gradient of the ambient parallel current and the gradient of the pressure (directly related to 
the diamagnetic current). The uncompensated perturbed magnetic energy is then set free 
in the narrow "inertial" (or "resonant" ) layer around r c and is converted to kinetic energy of 
poloidally rotating plasma (and drives reconnection when non-ideal effects are important). 
A kinetic description may be needed in order to address the physics of the inertial (resonant) 
layers. 

The internal kink instability has been known in plasma physics research from the begin- 
ning of the fusion programme. Already in 1973, it was considered as a plausible candidate 
to explain the q(0) > 1 tokamak stabilitiy criterion. It has been hypothesized that the 
internal kink mode can be involved in the disruptive instabilities. In Ref. jsj, an ideal Mag- 
netohydrodynamical (MHD) theory of the internal kink modes was developed. A nonlinear 
kinked neighboring equilibrium was found. However, the resulting nonlinear ideal-MHD 
amplitudes were too small to explain the disruptive phenomena in tokamaks although such 
qualitative properties as the negative voltage spikes and inward shifts in major radius pro- 
duced by the mode evolution agreed with the experimental observations. The next time the 
internal kink modes received increased attention was in 1975 when B. Kadomtsev proposed 
his model [3] for the sawtooth phenomena, discovered experimentally one year before |||. 
In this model, the growth of the magnetic island associated with the internal kink mode 
(poloidal mode number m = 1) plays a crucial role. Clearly, such a process is beyond the 
range of the ideal MHD theory considered in 

This has motivated the development of theoretical descriptions of the internal kink mode 
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which include various non-ideal effects. Thus, in Ref. js| the internal kink mode has been 
considered within a two-fluid model including resistivity, ion-ion collisions and diamagnetic 
effects (which were found to be stabilizing). Further development focused on a more detailed 
description of the resonant layer around the q — 1 magnetic surface. In Ref. [9|, the kinetic 
theory of the m = 1 internal instability was considered (kinetic electrons and fluid ions). 
Similarly to Ref. jg], a diamagnetic stabilization of the internal kink mode was found. In 
Ref. [lo| , the role of collisionless reconnection (driven by the electron inertia) on the evolution 
of the internal kink mode was studied in the regime 5 e ^> pi (here 5 e is the electron skin 
depth and pi is the ion gyroradius), neglecting ion Finite Larmor Radius (FLR) effects. Such 
a regime corresponds to very low values of plasma pressure (3 < m e /rrii. It has been pointed 
out, however, that an accurate treatment of the ion orbits (in terms of their gyro-average) 
is needed in order to consider realistic regimes with 5 e < pi. A technical difficulty arising 
here is due to the non-local character of the gyrokinetic polarization density which results 
in an integral quasineutrality equation (see also U[). This problem has been addressed in 
Ref. 12[ in the context of the tearing mode problem. The integral equation resulting from 
the non-local ion response has been solved by an iterative method (assuming the poloidal 
beta to be small). In Ref. 13|, an alternative approach to this problem has been suggested 
based on the Fourier expansion in the resonant layer with respect to the radial variable. This 
approach emphasizes the singular layer but retains the global nature of the mode through 
the boundary conditions. An interpolation formula has been used in order to treat the 
gyro-averages (it fails when the density or temperature gradients are too large). In 1991, 
Porcelli combined the approach of Ref. 13[ concerning the ions with electron inertia effects 
assuming, however, an adiabatic or isothermal electron temperature response (depending on 
the parameter regime). In Ref. [l^|, he found the so-called m = 1 collisionless tearing mode 
which dominates the ideal-MHD version of the internal kink mode when the ideal drive is 
sufficiently small (basically, when the ion gyroradius exceeds the ideal-MHD inertial length). 
This development has been used in the formulation of the revised sawtooth model 15[ by 
Porcelli et al in 1996 (this model takes into account the complex dynamics of the resonant 
layer which includes the ion FLR and reconnection effects). 

Despite the great progress achieved in the analytical understanding of internal kink 
modes, it is still an area of active research. Thus, very recently a unified theory of the 
internal kink and tearing modes has been developed in Ref. 16J. This theory provides an 
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accurate treatment of the ion orbits and the electron temperature response without assuming 
it to be adiabatic or isothermal. 

On the numerical side, there is a vast literature on simulations of internal kink modes 
using the fluid approach. For example, Ref. 17] analyzes the ideal stability of the inter- 



nal kink mode in shaped tokamak plasmas. Recent full-MHD simulations [18|, [19| studied 
various nonlinear regimes of the internal kink mode evolution and included some two-fluid 
(e. g. diamagnetic) effects. It is however desirable to have a first principles approach to the 
internal kink mode since the kinetic effects can be important in the resonant layer, and two- 
fluid models are often derived under a number of assumptions that may be too simplistic 



Such an attempt was undertaken for the first time with a gyrokinetic 20| particle-in-cell 21 1 



code in 1995. In Ref. 



22( , both the linear and nonlinear evolution of the internal kink mode 



were simulated, targeting the sawtooth collapse on a fast time scale. However, these early 
simulations operated at extremely small 8 < m e /mi and neglected ion FLR effects. This 
work has recently been continued in Ref. 23| (also at a very small plasma beta). In addition 
to particle-in-cell simulations, the internal kink mode has been considered by employing an 
eigenvalue approach to the solution of the gyrokinetic equations [24. |25|. 

In this paper, we study the linear evolution of the internal kink mode in straight-tokamak 
geometry employing the global gyrokinetic particle-in-cell method. Realistic values of plasma 
6 ^> m e /rrii are used. The objective of this paper is two-fold. First, we study how the kink 
mode properties change depending on the relation between the ion gyroradius and the MHD 
inertial length. Second, we explore, for the first time, what could be the limitations of the 
gyrokinetic PIC approach to this problem at realistic values of plasma beta. 

The structure of the paper is as follows. In Sec. [Til the basic equations and the discretiza- 
tion procedure are discussed. Simulations in the fluid and kinetic regimes are presented in 
Sec. IHIl Finally, our conclusions are summarized in Sec. IIVI 



II. BASIC EQUATIONS AND NUMERICAL APPROACH 



We use the linear two-dimensional 5f PIC-code GYGLES [26l428j|. The code is electro- 
magnetic and treats all particle species (ions and electrons) kinetically. It solves the gyroki- 
netic Vlasov- Maxwell system of equations (in the p» -formulation, see Ref. 20[ for details). 



The distribution function is split into a background part and a small time-dependent per- 
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turbation f s = F 0s + Sf s (the index s = i, e is used for the particle species). The background 
ion distribution function is taken to be a Maxwellian. The background electron distribution 
function is a shifted Maxwellian (to account for the equilibrium parallel current): 
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Here u = — n /(eno) and /i = (B/I?) ■ (V x B) with B the equilibrium magnetic field. 

Assuming the amplitude of the field perturbation to be small (this implies 5f s /F 0s <C 
1), the first-order perturbed distribution function can be found from the linearized Vlasov 
equation: 
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Here, [R/ \p||°' ) ] correspond to the unperturbed gyrocenter position and parallel velocity, and 
[R,W,pP] are the perturbations of the particle trajectories proportional to the electromag- 
netic field fluctuations [shown in Eqs. (jl])-([7j) below]. The perturbed part of the distribution 



function is discretized with markers (see Ref. 



29j for details): 
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where N p is the number of markers, 



(R v ,p\\ u , ijl v ) are the marker phase space coordinates 



and w su is the weight of a marker. The equations of motion are 
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with <fi and An being the perturbed electrostatic and magnetic potentials, \x the magnetic 
moment, m the mass of the particle, = b • V x A*, b* = V x A* /Bp A* = A + 
(mp\\/q)b the so-called modified vector potential, A the magnetic potential corresponding 
to the equilibrium magnetic field B = V x A and b = J5/B the unit vector in the direction 
of the equilibrium magnetic field. The gyro-averaged potentials are defined as usual: 



(0) = /-0(R + p), (A ll ) = f-A^R + p) 



(8) 



where p is the gyroradius of the particle and 9 is the gyro-phase. Numerically, the gyro- 
averages are computed sampling a sufficient number of the gyro-points on the gyro-ring 
around the gyrocenter position of the marker 
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30]. 



The perturbed electrostatic and magnetic potentials are found self-consistently from the 
gyrokinetic quasineutrality equation and parallel Ampere's law 3l| : 

QiFoi , 
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where 5n s = J d 6 Z Sf s 5(H + p — x) is the gyrocenter density, 6j\\ a = q s J d 6 Z 5f s p\\ 5(R + 
p — x) is the gyrocenter current, q s is the charge of the particle, d 6 Z = B% dRdpy d/id^ is 
the phase-space volume, p s = y/m s T s / (eB) is the thermal gyroradius and f3 s = ^qtiqT s /BI is 
the plasma beta corresponding to a particular species. Note that the polarization density in 
Eq. (Q is given by an integral operator and includes a non-local effect of the ion gyro-orbit. 

We have found, however, that the computations can be simplified by replacing the 
quasineutrality condition, Eq. with 



fc(r) 
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with r c being the position of the singular layer. This representation implies that the exact 
expression for the polarization density is used around the resonant position q(r c ) = 1 (where 
the relevant radial scale can be smaller than the ion gyroradius), whereas the long- wavelength 
approximation is used outside the resonant layer (which is justified since the kink mode 
structure is global in the ideal region). The effect of this representation is to replace (in 
the ideal region only!) the kinetic Alfven waves by shear Alfven waves, which have more 
favourable (for numerics) properties experiencing more physical damping on small radial 
scales. In our simulations, we have used Ar c = 0.2r a and p c = 8. For this choice of the 
parameters, pi Ar c in all cases considered. 

The electrostatic and magnetic potentials are discretized with the finite-element method 
(Ritz-Galerkin scheme): 

N a N s 

0(x) = ]T Mz(x) , A„ (x) = £ aiAi(x) , (13) 
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where A;(x) are finite elements (tensor products of B splines |32|, [33[ ) , N s is the total number 
of the finite elements, <pi and a/ are the spline coefficients. The numerical treatment of the 
nonlocal gyrokinetic polarization density has been described in Ref. 30(. Homogeneous 
Dirichlet boundary conditions are applied for <j) and An both on the axis and at the plasma 



edge. Further details of the numerical discretization can be found in Refs. 



29 



34 



III. SIMULATIONS 

We consider a straight tokamak (a screw pinch) with minor radius r a and "major radius" 
Ro (the pinch is a topological torus with the length L = 2ttRq and periodic boundary 
conditions along the axis). The plasma consists of hydrogen ions and electrons (with a 
realistic mass ratio). The safety factor is given by the expression q(r) = qo + (1 — qo)(r/r c ) Pq 
(with q = 0.6 and p q = 1 in our simulations). The background magnetic field is determined 
by the MHD pressure balance condition: 

r(* + ^) + - = °> b -tV b - < 14 > 

dr \ 2fx J /i r q{r)Ro 

This equation is solved with the "initial condition" B z (r = 0) = Bq. It can be seen that 
B z (r) « B often gives a good approximation for the equilibrium magnetic field (when the 
effect of the plasma pressure is small). 

The plasma temperature and density profiles are given by the expressions: 



77, (s) = N Q exp 
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(15) 



where s = r/r a with r a the radius of the pinch. The shape of these profiles can be flexibly 
tailored by adjusting the parameters At and A n (which determine the profile width), st 
and s n (position of the maximal gradient), Lt s and L n (gradient lengths). 

The numerical resolution needed in the simulations depends on the physical parameters. 
It has been observed that the radial resolution must always be sufficient to resolve the ion 
gyroradius (it is required since the Alfven continuum must be resolved). For cases when 
the electron inertia is of importance, the electron skin depth must also be resolved. In the 
"poloidal" direction, N$ = 4 splines have been found sufficient since only one poloidal mode 
is kept in the linear straight-tokamak simulations. The marker resolution is kept on the level 
of 700 — 1000 markers per grid cell. 



A. Internal kink mode: "MHD regime" 



In this subsection, we consider the internal kink instability in the MHD regime. Physi- 
cally, it implies that the width of the resonant layer (which in this case coincides with the 
MHD inertial- layer width A#) is much larger than the ion gyroradius pi. The inertial-layer 
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width is given by the expression 
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Here, q' = dq/dr, k z = n/ Rq and k ■ B = B z (m — qri)/(qRo). For the internal kink mode, 
the poloidal mode number m = 1 and the toroidal mode number n = 1. 

Analytical ideal-MHD theory gives the following expression for the growth rate: 

ir a = q'{r c )\ H , t a = R /v A . (18) 

Here, v A = B z j JpZ^m^iQ is the Alfven velocity. Another way (used here) is to solve numer- 
ically (with the shooting method) the ideal-MHD eigenvalue problem: 

A (j/io^noT 2 + (k • B) 2 ]r 3 |^ - g x {r)i = (19) 

employing the boundary conditions £'(r = 0) = and £(r = r a ) = 0. Here, ^(r) is the MHD 
displacement and £' = d^/dr. 

We consider a straight tokamak with "major radius" Rq = 5 m, minor radius r a — 1 m, 
magnetic field Bq = 2.5 T, plasma temperature T« = T e = 5 keV, and plasma density 
Nq = 2 x 10 19 m? (which corresponds to /? = 0.0128 at the temperature chosen). Both the 
plasma density and the temperature are flat. For these parameters, an unstable internal 
kink mode exists (destabilized by the gradient of the ambient parallel current). 

We consider a sequence of straight-tokamak equilibria corresponding to different locations 
of the rational flux surface r c . In Fig. [H the ideal-MHD inertial-layer width is shown as 
a function of r c . One sees that A# exceeds both the ion gyroradius pi = 2.8 x 10~ 3 m and 
the electron skin depth 5 e = p e /\fWe — 1-7 x 10~ 3 m for all equilibria considered. 
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In Fig. [21 the internal kink mode growth rate resulting from the gyrokinetic PIC sim- 
ulations is compared to the results of the ideal-MHD eigenvalue calculations [the MHD 
eigenvalue problem Eq. ( I19p has been numerically solved using the shooting method]. One 
sees that the mode is more unstable for larger r c , i. e. when a larger plasma column is 
involved in the instability [this can also be seen formally from Eq. (Il6p ]. The quantitative 
agreement between the ideal MHD and the gyrokinetic simulations is very good when the 
ion temperature is small Tj = 200 keV (implying ion gyroradius small). When the ion 
gyroradius increases, the kink mode becomes less unstable (although the quialitative de- 
pendence on r c remains the same). The mode is still an MHD-like internal kink instability 
but this will change when k±pi > 1 with k± ~ 1/A# (see Sec. IIII Bp . One can quantify the 
ion-FLR effects considering a sequence of internal kink modes keeping all the parameters 
constant except the plasma temperature (Tj = T e ). Note that the ideal-MHD growth rate 
does not depend on the temperature (for a flat profile such as used in these simulations). In 
contrast, the gyrokinetic simulations show that the kink mode growth rate decreases with 
the plasma temperature (see Fig. [3]). This FLR-stabilization effect can be quite substantial 
when the ion temperature is large enough. Note that a similar ion-FLR stabilization has 
already been observed for other MHD modes (see the gyrokinetic simulations of the Toroidal 
Alfven Eigenmodes in Ref. 27|). The stabilization is caused by the gyro-averaging operation 
acting on the perturbed electromagnetic field. Clearly, this effect is absent in the ideal MHD 
description and can be found only with a kinetic treatment. 

In Fig.Hl we compare the ideal-MHD eigenmode structure with the radial pattern found in 
the gyrokinetic simulations. An excellent agreement is found. The ideal-MHD displacement 
£ ~ 0/r corresponding to Fig. [4] has the well-known top-hat structure. The poloidal fluid 
velocity vg ~ d<p/dr is clearly strongly increased at the rational flux surface, as must be the 
case for the internal kink mode. 

In Fig. [U the internal kink mode evolution is shown (here, the rational flux surface 
is located at r c = 0.7; the plasma temperature Tj = T e = 5 keV). One sees how the 
initial perturbation (which has had a Gaussian shape in the radial direction) decays in the 
continuum of the shear Alfven waves and is then re-organized as an unstable internal kink 
mode. The ion gyroradius must be resolved in order to correctly reproduce the continuum- 
decay phase of the mode evolution. The frequency of this mode is zero (this is a well-known 
property of the ideal MHD modes, related to the Hermitian symmetry of the underlying 
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equations; this property can also be reproduced in the gyrokinetic simulations). 

So far, we have considered the current- driven internal kink mode (the plasma pressure 
was chosen to be flat). Now, let us include the pressure destabilization in our simulations. 
In Fig. [6J the kink- mode growth rate is plotted as a function of the plasma density gradient. 
For the density profile, we choose s n = 0.5, A n = 0.2 [see Eq. (flol) ]. The plasma temperature 
Tj = T e = 3 keV is taken to be flat. The resonant flux surface has been located at r c = 0.7. 
Two cases, one with a large magnetic field B = 5 T and the other with a moderate magnetic 
field B = 2.5 T, were considered [keeping /3(s n ) = 0.00773 constant]. One sees from Eq. (fI8|) 
that the ideal-MHD result does not depend on the magnetic field strength. Indeed, the 
gyrokinetic simulations reproduce this property. The agreement between the gyrokinetic 
simulations and the MHD result is good for both values of the ambient magnetic field. The 
kink- mode growth rate increases with the density (pressure) gradient, as expected. 

Summarizing, we have considered the internal kink mode in the regime pi < Xh (ions in 
the resonant layer are magnetized). In this regime, the properties of the kink mode agree 
well with the ideal MHD expectations (thus providing a benchmark for our simulations). 
The only effect found beyond ideal MHD, is an ion-FLR stabilization of the internal kink 
mode (which however can be quite substantial). 



B. Internal kink mode: "FLR regime" 

We now consider the case when the thermal ion gyroradius exceeds the ideal-MHD inertial 
width pi ^> A#. The parameters here are as follows. The straight tokamak has the minor 
radius a = 0.5 m, the ambient magnetic field Bq = 1 T, and the plasma density uq = 
1.6 x 10 18 m~ 3 (it corresponds to /? = 0.00644 when Tj = T e = 5 keV). Equilibria with the 
major radii Rq = 5 m and Rq = 20 m have been considered (changing the major radius 
one can change the aspect ratio and, consequently, the inertial-layer width Xh while pi/r a is 
kept fixed). The comparison between the inertial layer width Xh, the thermal ion gyroradius 
Pi, and the electron skin depth 5 e is shown for different aspect ratios in Fig. [7J One sees 
that the ion gyroradius is indeed much larger than all other relevant scales, especially when 
Ro = 20 m (this case corresponds to a particularly weak, almost vanishing ideal-MHD drive). 



In this regime (ions are demagnetized in the resonant layer), Ref. 14( predicts the exis- 



tence of an unstable collisionless m — 1 tearing mode. In the case Xh — > 0, the growth rate 
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of this mode is given by the expression 



s,(r c )v A (r c ) (Serf) 1 ' 3 Bl . r dq 

J = —, — • VA= VJ^' s, = 5^ (20) 

In Fig. [HI we plot the growth rate resulting from the gyrokinetic PIC simulations (correspond- 
ing to the pinch with Rq = 20 m, i. e. very small A#) compared to the analytic prediction 
for the collisionless m — 1 tearing mode as a function of the resonant layer position. In 
contrast to the "ideal-MHD" regime (see Fig. [2]), the mode is more unstable for smaller r c , 
indicating a change in the underlying physical mechanism: reconnection (electron physics) 
in addition to the poloidal plasma rotation (ions). The agreement between the gyrokinetic 
result and the theoretical prediction Eq. (120]) is very good. 

In Fig. [HI we consider the dependence of the kink-mode growth rate on the ion tem- 
perature (ion-FLR effect). One sees that in the regime of demagnetized ions (inside the 
resonant layer), the mode is further destabilized when the ion gyroradius increases. This 
FLR-destabilization is in contrast to the MHD-type internal kink mode considered previ- 
ously, which was FLR-stabilized (see Fig. [3]). The numerical result has been compared with 
the analytical expression Eq. ( 1201) . The agreement is again very good. For comparison, we 
plot the growth rate predicted by the ideal MHD theory, which is more than two orders of 
magnitude smaller. It indicates that the kinetic effects can dominate the kink mode physics, 
certainly around ideal marginal stability. 

Finally, we consider the effect of the electron temperature gradient on the collisionless 
m = 1 tearing mode. We employ a temperature profile with St = r c /r a = 0.5 (the position 
of the maximal gradient coincides with that of the resonant flux surface) and A s = 0.2 [see 
Eq. (115j) ]. The ion temperature and the plasma density are taken to be flat. As a consequence 
of the diamagnetic effect associated with the finite electron temperature gradient, the mode 
acquires a finite frequency and becomes a drift-tearing instability (see e. g. Q,Q). The 
frequency resulting from the gyrokinetic PIC simulations is plotted in Fig. [10] as a function 
of Lrpe (electron temperature gradient length) for two different values of the major radius: 
Rq = 5 m and Rq = 10 m (recall that the ideal-MHD drive and, consequently, Xh decrease 
with the aspect ratio). One sees that the frequency is the same for both values of Rq. Hence, 
it is set exclusively by the diamagnetic frequency at the position of the resonant flux surface 
WTe = l/(-B r c) dT e /dr. In fact, the fre quen cy of the collisionless m — 1 drift-tearing mode 



appears to satisfy uj ~ 0.4 u Te (cf. Ref. jl2|). 
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The growth rate as a func 



;ion of Lt e (electron temperature gradient length) is plotted in 



Fig. [TTJ As expected [12|, |l6| , the electron temperature gradient stabilizes the collisionless 



drift-tearing m = 1 mode. Physically, it has been shown in Ref. 16| that the drift-tearing 



mode can couple to a stable Kinetic Alfven Wave (KAW) which leads to a stabilization of 



the mode. Another way to explain this stabilization effect has been elaborated in Ref. 3?| 
(in the fishbone context). In short, the finite-frequency mode can interact with the Alfven 



continuum (or, alternatively, with the KAWs as in [16[) and thus undergo continuum damp- 
ing (as is well-known from the context of Toroidal Alfven Eigenmodes 38M40I] ) . Of course, 
the frequency of the drift-tearing m — 1 mode is quite small, so that the continuum damp- 
ing acts only very close to the resonant flux surface (the surface of vanishing k\\) where the 
condition u 2 = k 2 (r a) v a{ t a) can be satisfied. The frequency of the m — 1 drift-tearing 
mode increases with the electron temperature gradient which makes the continuum damp- 
ing more efficient at smaller L? e - to the point of a complete stabilization. Note that in 
Fig. [TT] the stabilization appears to be stronger for the case with the larger major radius 
Rq = 10 m. In this case, the absolute value of the growth rate would be smaller compared 
with the configuration with Rq = 5 m, even for the flat temperature profiles, because the 
growth rate inversely scales with the Alfven time ta ~ Ro- Thus, the continuum damping, 
which is determined by the mode frequency (set by the electron temperature gradient and 
thus equal for both values of Rq), corresponds to a larger fraction of the smaller growth rate 
when Rq = 10 m. 

In Fig. [121 the radial mode structure is shown. One sees that in addition to the conven- 
tional internal kink eigenmode, a complicated fine-scale structure appears at the resonant 
flux surface. This structure is caused by the continuum damping of the instability. One can 
see that it becomes more pronounced when the electron temperature gradient increases (see 
Fig. [13]) . We can estimate the position of the shear Alfven resonances using the expression: 

r A = r c ± , t a = Rq/va (21) 

q'[r c ) 

For the safety factor profile chosen, q'(r c ) = 0.8. Fig. [13] indicates that the fine-scale structure 
developed could be understood as Kinetic Alfven Waves (recall that we use a non-local 
expression for the polarization density) resonantly excited by the collisionless m = 1 tearing 
mode in the positions approximately satisfying the shear Alfven wave resonance condition 
uj 2 ~ k 2 (r A )v A (rA) (the Kinetic Alfven Wave resonance would be more precise). The drift- 
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tearing (kinetic-kink) mode is damped because of its coupling to the Kinetic Alfven Waves 
(a "continuum" ) at the resonant positions. A finite frequency of the tearing mode is required 
for this process to function. This frequency is provided by the diamagnetic effect due to a 
finite electron-temperature gradient. 



IV. CONCLUSION 



In this paper we have studied the internal kink modes in a straight-tokamak geometry 
using the global gyrokinetic particle- in-cell code GYGLES. Both electron and ion gyrocenters 
were treated kinetically, but collisions were ignored. The simulations have shown that the 
kink mode properties depend strongly on the ratio between the ideal-MHD inertial-layer 
scale and the gyroradius. In the "MHD regime", the kink mode becomes more unstable 
if the rational flux surface "moves" (during the parameter scan) outwards. In the "FLR 
regime", however, the kinetic- kink mode (the "collisionless m = 1 tearing mode" Jjj]) 
is more unstable for resonant magnetic surfaces located closer to the magnetic axis (l/r c 
dependence of the growth rate). Similarly, the scaling with respect to the ion temperature 
is also opposite: the FLR-stabilization for the MHD-type kink mode has been observed 
whereas the "collisionless m = 1 tearing mode" is FLR- destabilized. One could speculate 
[ijj] that the most unstable m — l,n — 1 instabilities observed in real hot plasmas often 
correspond to the "collisionless m — 1 tearing" (kinetic-kink) modes rather than the classical 
ideal-MHD internal kink modes since usually the real plasmas are at the marginal MHD-kink 
stability boundary so that the ion gyroradius can easily overcome the ideal inertial length 
(thus bringing the instability into the kinetic regime). The kinetic- kink modes, however, 
can be sensitive to the shape of plasma profiles and can be stabilized e. g. by the electron 
temperature gradient (through a combination of the diamagnetic and continuum-damping 
effects) . 

Looking forward, gyrokinetic simulations of the internal kink modes in tokamak geometry 
would be of interest. In these simulations, the effects of the guiding center orbits in the 
resonant region, the kinetic trapped-ion and fast-particle (e. g. fishbone) effects in the ideal 
region should be addressed. Also, gyrokinetic simulations of interchange instabilities (both 
in straight and toroidal geometries) could be performed in the future. 
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FIG. 1: (Color online) Ideal-MHD inertial-layer width Ah as a function of r c [with k\\(r c ) = 0] 
compared to the ion thermal gyroradius pi (computed for = 5 keV) and to the electron skin 
depth 5 e . Screw pinch geometry with i?o = 5 m, a = 1 m, Bq = 2.5 T is considered. For these 
parameters, the fluid inertial length is larger than the kinetic radial scales p% = 2.8 x 10 -3 m and 
5 e = 1.7 x 1(T 3 m. 
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FIG. 2: (Color online) Internal kink mode growth rate resulting from the gyrokinetic PIC simu- 
lations (at two different ion temperatures) compared to the results of the ideal-MHD eigenvalue 
calculations (the MHD eigenvalue problem has been numerically solved using the shooting method). 
The agreement between the ideal MHD and the gyrokinetic simulations is very good, especially 
at the smaller ion temperature (i. e. at the smaller ion gyroradius). One sees that the gyrokinetic 
kink mode is somewhat stabilised when the ion gyroradius increases. This FLR effect is absent in 
the ideal MHD description. 
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FIG. 3: (Color online) Internal kink mode growth rate (resulting from the gyrokinetic PIC simu- 
lations) as a function of the plasma temperature (Tj = T e ). The growth rate decreases at larger Tj 
which should be attributed to FLR-stabilization. The rational flux surface is localized at r c /a = 0.5. 
For comparison, the ideal-MHD result is shown (which does not depend on the ion temperature). 
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FIG. 4: (Color online) Internal kink mode radial structure in the regime with A# 3> p^. The 
mode changes abruptly at the rational flux sruface r c = 0.5 which is a well-known property of the 
intenal kink modes. One sees a perfect agreement between the eigenmode structures obtained from 
the initial-value gyrokinetic PIC simulations and solving (numerically) an ideal-MHD eigevalue 
problem. 
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FIG. 5: (Color online) Internal kink mode evolution in the regime with Xh 3> p%- The rational flux 
sruface is located at r c = 0.7. The plasma temperature = T e = 5 keV. The initial perturbation 
(a Gaussian) decays in the continuum of the shear Alfven waves which are reorganized afterwards 
into the internal kink eigenmode. 
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FIG. 6: (Color online) Pressure gradient effect on the kink mode in the MHD regime. Equilibria 
with "large" and "moderate" magnetic fields (5 = 5T and B = 2.5 T) have been compared. The 
plasma temperature T, = T e = 3 keV. The result does not depend on B (as expected according to 
MHD theory) and the agreement between the gyrokinetic simulations and the MHD computation 
is good. 
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FIG. 7: (Color online) Ideal-MHD inertial-layer width compared with the ion thermal gyroradius 
and the electron skin depth. Screw pinches with a = 0.5 m, Bq = 1 T, Rq = 5 m and Rq = 20 m are 
considered. Here, the MHD drive is much smaller compared to the case shown in Fig. [T] (especially 
when Rq = 20 m). As a consequence, the ion gyroradius is much larger than the MHD inertial 
length indicating the importance of the FLR effects (and ion sub-Larmor scales) in the mode 
evolution (kinetic regime of the m = 1, n = 1 mode). 
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FIG. 8: (Color online) Internal kink mode growth rate resulting from the gyrokinetic PIC sim- 
ulations compared to the analytic prediction for the collisionless m = 1 tearing mode [14J] 
7 ~ ((5 e p?) 1//3 /r c . Here, the parameters chosen are Tj = T e = 11 keV and Rq = 20 m. One 
sees a very good agreement between the theory and simulations. 
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FIG. 9: (Color online) Internal kink mode growth rate (resulting from the gyrokinetic PIC simula- 
tions for Rq = 20 m) as a function of the plasma temperature (Tj = T e ). The growth rate increases 
with the temperature contrary to the ideal-MHD case (see Fig. [3]). The rational flux surface is 
localized at r c = 0.4. The agreement of the simulations with the theory is very good. For reference, 
the ideal-MHD growth rate corresponding to the parameters chosen is plotted. One sees that the 
growth rate of the kinetic-kink (reconnecting) mode is more than two orders of magnitude larger 
than the growth rate of the marginally-unstable ideal-MHD mode. 
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FIG. 10: (Color online) Frequency of the kinetic internal kink mode (collisionless m = 1 tearing 
mode) as a function of the electron temperature gradient. Ion temperature and density are kept 
constant. For comparison, an estimate for the drift-tearing mode frequency u = 0.4 wx e [3] is 
plotted. One sees that the agreement is good. The mode properties for two different aspect ratios 
(and, consequently, \h) are compared (major radius is varied whereas the minor radius is kept 
fixed). One sees that the mode frequency does not depend on Xh and, hence, is set exclusively by 
the diamagnetic frequency associated with the electron temperature gradient. 
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FIG. 11: (Color online) Growth rate of the kinetic internal kink mode (collisionless m = 1 tearing 
mode) as a function of the electron temperature gradient (same parameters as in Fig. [TO]) . In 
contrast to the frequency, the growth rate is much smaller at larger Rq (note the normalization 
to Tj\ oc 1 / Rq . ) In this case, the continuum damping (set by the mode frequency which is equal 
for both values of i?o) corresponds to a larger fraction of the smaller growth rate, making the 
stabilization mechanism more effective. 
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FIG. 12: (Color online) Kinetic-kink eigenmode (Rq = 10 m) compared at different electron tem- 
perature gradient lengths Lxe- A fine-scale structure develops around the rational flux surface 
[indicating the continuum damping at the position of the resonance oj = ±k\\{r a)v a, see Fig. [13] 
for details]. The fine-scale structure is more pronounced (see Fig. [T3|) at larger electron tempera- 
ture gradients (where the mode frequency is larger and, consequently, the continuum damping is 
stronger). 
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FIG. 13: (Color online) The eigenmodes from Fig. [12] zoomed around the rational flux surface. 
For comparison, the ion gyroradius and electron skin depth scales are plotted. Also, the positions 
of the shear Alfven resonances are indicated (the kinetic Alfven wave resonances would be more 
precise). One sees that the relevant radial scales are smaller than the ion gyroradius. Kinetic 
Alfven waves are excited at the resonant positions causing a stabilization of the collisionless m = 1 
drift-tearing mode (through the "continuum" damping mechanism). 
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